#include <boost/multiprecision/cpp_int.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>
#include <iostream>
#include <fstream>
#include <vector>
#include <algorithm>
#include <cmath>
#include <cstdlib>

struct line{
	boost::multiprecision::cpp_rational a,b,c;
	int x1,x2,y1,y2;
	
	line(boost::multiprecision::cpp_rational a, boost::multiprecision::cpp_rational b, boost::multiprecision::cpp_rational c, int x1, int x2, int y1, int y2): a(a), b(b), c(c), x1(x1), x2(x2), y1(y1), y2(y2) {}
		
	friend bool operator <(const line &L1, const line &L2);
};

void Perwitt(cv::InputArray src, cv::OutputArray dst, double alpha = 1.0, double beta = 1.0, double gamma = 0.0);
double get_real_area(std::vector<std::vector<cv::Point> > contours, int idx, cv::Size size);
void split(std::vector<std::vector<cv::Point> > contours, int idx);
void draw_line(const line &L, cv::InputOutputArray img, int RGB_R = 255, int RGB_G = 255, int RGB_B = 255);

const double PI = acos(-1);
const double min_length = 50;
const double min_area = 100;
const double min_length2 = 10;
const double min_delta_theta = PI * 65 / 180;
const double ang_eps = PI / 18;

std::ofstream logfile;
std::vector<line> lines;
double Rgh = 0, RghWgt = 0, RghMax = 0, RghMin = 1e8;
bool logfile_enabled, output_enabled;

int main(int argc, char *argv[])
{
	std::ios::sync_with_stdio(false);
	if (argc < 2 || argc > 4)
	{
		std::cerr<< "Usage: roughness input_file [output_file] [log_file]\n\n";
		return 1;
	}
	output_enabled = (argc >= 3);
	logfile_enabled = (argc >= 4);
	srand(30970); /* 30970 is a randomly selected seed */
	if (logfile_enabled) logfile.open(argv[3], std::ios::app);
	cv::Mat input = cv::imread(argv[1]), tmp01, tmp02, tmp03, tmp04, tmp05, tmp06 = cv::Mat::zeros(input.size(), CV_8UC3);
	cv::Mat tmp07 = cv::Mat::zeros(input.size(), CV_8UC3);
	cv::GaussianBlur(input, tmp01, cv::Size(3, 3), 0);
#ifdef DEBUG
	cv::imwrite("./tmp/tmp01.png", tmp01);
#endif
	Perwitt(tmp01, tmp02);
#ifdef DEBUG
	cv::imwrite("./tmp/tmp02.png", tmp02);
#endif
	cv::cvtColor(tmp02, tmp03, cv::COLOR_RGB2GRAY);
#ifdef DEBUG
	cv::imwrite("./tmp/tmp03.png", tmp03);
#endif
	cv::addWeighted(tmp03, 0.0, tmp03, 2.0, -20, tmp04);
#ifdef DEBUG
	cv::imwrite("./tmp/tmp04.png", tmp04);
#endif
	cv::threshold(tmp04, tmp05, 0, 255, cv::THRESH_OTSU);
#ifdef DEBUG
	cv::imwrite("./tmp/tmp05.png", tmp05);
#endif
	std::vector<std::vector<cv::Point> > contours;
	std::vector<cv::Vec4i> hierarchy;
	cv::findContours(tmp05, contours, hierarchy, cv::RETR_LIST, cv::CHAIN_APPROX_SIMPLE);
	for (int i = 0; i >= 0; i = hierarchy[i][0])
		if (cv::arcLength(contours[i], true) >= min_length && get_real_area(contours, i, input.size()) >= min_area)
		{
			cv::drawContours(tmp06, contours, i, cv::Scalar(255, 255, 255));
			split(contours, i);
		}
	if (logfile_enabled)
	{
		logfile << "ROUGHNESS(avg. ) = " << Rgh / RghWgt << "\n";
		logfile << "ROUGHNESS(max. ) = " << RghMax << "\n";
		logfile << "ROUGHNESS(min. ) = " << RghMin << "\n";
	}
	std::cout << Rgh / RghWgt << "\n";
#ifdef DEBUG
	cv::imwrite("./tmp/tmp06.png", tmp06);
#endif
	for (auto v: lines) draw_line(v, tmp07, rand() & 255, rand() & 255, rand() & 255);
	if (output_enabled) cv::imwrite(argv[2], tmp07);
	if (logfile_enabled) logfile.close();
	return 0;
}

void Perwitt(cv::InputArray src, cv::OutputArray dst, double alpha, double beta, double gamma)
{
	cv::Mat PerwittL = (cv::Mat_<double>(3,3) << -1.0, 0.0, 1.0, -1.0, 0.0, 1.0, -1.0, 0.0, 1.0), 
	PerwittR = (cv::Mat_<double>(3,3) << 1.0, 0.0, -1.0, 1.0, 0.0, -1.0, 1.0, 0.0, -1.0), 
	PerwittU = (cv::Mat_<double>(3,3) << -1.0, -1.0, -1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0), 
	PerwittD = (cv::Mat_<double>(3,3) << 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, -1.0, -1.0, -1.0);
	cv::Mat tmp01, tmp02, tmp03, tmp04, tmp05, tmp06;
	cv::filter2D(src, tmp01, -1, PerwittL);
	cv::filter2D(src, tmp02, -1, PerwittR);
	cv::filter2D(src, tmp03, -1, PerwittU);
	cv::filter2D(src, tmp04, -1, PerwittD);
	cv::addWeighted(tmp01, 1.0, tmp02, 1.0, 0, tmp05);
	cv::addWeighted(tmp03, 1.0, tmp04, 1.0, 0, tmp06);
	cv::addWeighted(tmp05, alpha, tmp06, beta, gamma, dst);
}

double get_real_area(std::vector<std::vector<cv::Point> > contours, int idx, cv::Size size)
{
	cv::Mat tmp = cv::Mat::zeros(size, CV_8UC1);
	cv::drawContours(tmp, contours, idx, cv::Scalar(255, 255, 255), cv::FILLED);
	return cv::countNonZero(tmp);
}

inline int __sqr(int x)
{
	return x*x;
}

void __LSM(std::vector<cv::Point> points, int begin, int end)
{
	int n = end - begin + 1, x1 = points[begin].x, x2 = points[begin].x, y1 = points[begin].y, y2 = points[begin].y;
	double length = 0;
	boost::multiprecision::cpp_int sum_x = 0, sum_y = 0, sum_xy = 0, sum_xx = 0, sum_yy = 0;
	for (int i = begin; i <= end; i++)
	{
		sum_x += points[i].x;
		sum_y += points[i].y;
		sum_xy += points[i].x * points[i].y;
		sum_xx += points[i].x * points[i].x;
		sum_yy += points[i].y * points[i].y;
		x1 = std::min(x1, points[i].x);
		x2 = std::max(x2, points[i].x);
		y1 = std::min(y1, points[i].y);
		y2 = std::max(y2, points[i].y);
		if (i < end) length += sqrt(__sqr(points[i + 1].x - points[i].x) + __sqr(points[i + 1].y - points[i].y));
	}
	boost::multiprecision::cpp_rational a = sum_xy - boost::multiprecision::cpp_rational(sum_x * sum_y, n), b = sum_xx - boost::multiprecision::cpp_rational(sum_x * sum_x, n), c, sigma2, delta;
	if (b == 0)
	{
		b = a;
		a = sum_yy - boost::multiprecision::cpp_rational(sum_y * sum_y, n);
	}
	if (a == 0 && b == 0) return;
	c = - a * boost::multiprecision::cpp_rational(sum_x, n) - b * boost::multiprecision::cpp_rational(sum_y, n);
	for (int i = begin; i <= end; i++)
	{
		delta = a * points[i].x + b * points[i].y + c;
		sigma2 += delta * delta / (a * a + b * b);
		if (double(abs(delta)) / sqrt(double(a * a + b * b)) > RghMax) RghMax = double(abs(delta)) / sqrt(double(a * a + b * b));
		if (double(abs(delta)) / sqrt(double(a * a + b * b)) < RghMin) RghMin = double(abs(delta)) / sqrt(double(a * a + b * b));
	}
	sigma2 /= n - 1;
	if (logfile_enabled) logfile << a << "x+" << b << "y+" << c << "=0\tsigma^2=" << sigma2 << "\n";
	lines.push_back(line(a, b, c, x1, x2, y1, y2));
	Rgh += sqrt(double(sigma2)) * length;
	RghWgt += length;
}

inline double __get_angle(cv::Point p1, cv::Point p2)
{
	return (p1.x == p2.x)?(PI / 2):atan((double)(p2.y - p1.y) / (p2.x - p1.x));
}

void split(std::vector<std::vector<cv::Point> > contours, int idx)
{
	int n = contours[idx].size();
	double ang1, ang2, theta;
	std::vector<cv::Point> points;
	std::vector<int> list;
	for (int i = 0; i < n; i++) points.push_back(contours[idx][i]);
	for (int i = 0; i < n; i++) points.push_back(contours[idx][i]);
	n *= 2;
	list.push_back(0);
	
	double len1, len2;
	int pre, nxt;
	for (int i = 1; i < n - 1; i++)
	{
		len1 = len2 = 0;
		pre = nxt = i;
		while (pre > 0 && len1 < min_length2)
		{
			pre--;
			len1 += sqrt(__sqr(points[pre].x - points[pre + 1].x) + __sqr(points[pre].y - points[pre + 1].y));
		}
		while (nxt < n - 1 && len2 < min_length2)
		{
			nxt++;
			len2 += sqrt(__sqr(points[nxt].x - points[nxt - 1].x) + __sqr(points[nxt].y - points[nxt - 1].y));
		}
		ang1 = __get_angle(points[i], points[pre]);
		ang2 = __get_angle(points[i], points[nxt]);
		theta = fabs(ang1 - ang2);
		if (theta > PI / 2) theta = PI - theta;
		if (theta > min_delta_theta)
		{
			list.push_back(i);
			if (i >= n / 2) break;
		}
	}
	if (list[list.size() - 1] < n / 2) list.push_back(n - 1);
	
	n = list.size();
	for (int i = 0; i < n - 1; i++) __LSM(points, list[i], list[i + 1]);
}

inline void draw_line(const line &L, cv::InputOutputArray img, int RGB_R, int RGB_G, int RGB_B)
{
	if (L.a == 0 && L.b == 0) return;
	if (abs(L.a) < abs(L.b)) cv::line(img, cv::Point(L.x1, int((- L.a * L.x1 - L.c) / L.b)), cv::Point(L.x2, int((- L.a * L.x2 - L.c) / L.b)), cv::Scalar(RGB_B, RGB_G, RGB_R));
	else cv::line(img, cv::Point(int((- L.b * L.y1 - L.c) / L.a), L.y1), cv::Point(int((- L.b * L.y2 - L.c) / L.a), L.y2), cv::Scalar(RGB_B, RGB_G, RGB_R));
}
